clear all
tempfile tempsave
set seed ${seed}

do "${scripts}_output.do"
do "${scripts}_texgraph.do"

/*Macros*/
local data 		"National Longitudinal Study of Adolescent to Adult Health"
local source 	"Own calculations."
local signif 	"Significance levels:" ///
				"* \(p<0.10\), ** \(p<0.05\), *** \(p<0.01\)."
local error 	"Standard errors are clustered at the school level."
local controls 	"\emph{Child Controls}:" ///
				"Firstborn dummy, linear birth cohort trend (in months) by gender," ///
				"20 principal components of the full matrix of genetic data." ///
				"\emph{Family Controls}:" ///
				"Age of mother at birth, years of education of both mother and father," ///
				"average potential wages of both mother and father," ///
				"the standard deviation of potential wages of both mother and father, dummies for non-US born mothers and fathers," ///
				"a dummy for Christian religion, state fixed effects."
local controlfct	"\emph{Control Function}: Share white, share single mothers, maternal education (average),"  ///
					"share females, share migrants. All control function variables are calculated"  ///
					"as leave-cohort-out school averages."
local standard  "All right-hand side variables are standardized on the estimation sample ($\mu=0$, $\sigma=1$)."
local notesize scriptsize

local format pdf
set scheme white_tableau
graph set window fontface "Palatino Linotype"
colorpalette viridis, nograph
local c1="`r(p1)'"
local c2="`r(p5)'"
local c3="`r(p10)'"
local c4="`r(p15)'"

// -----------------------------------------------------
// Figure 1
// -----------------------------------------------------
local name fig1
use "${data}data_estim.dta" if sample1==1, clear

qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

foreach x in pgs_edu f_over1{
	if "`x'"=="pgs_edu"{
		local gtitle "PGI{superscript:EA}"
		local title_nc "No Controls"
		local title_wc "Incl. Controls"
	}
	if "`x'"=="f_over1"{
		local gtitle "Q{superscript:S}"
		local title_nc " "
		local title_wc " "
	}

	/*Uncontrolled*/
	reghdfe eduyears `x' [aw=w0], absorb(w1state) cluster(HsId)
	local b=string(_b[`x'],"%9.3f")
	local se=string(_se[`x'],"%9.3f")
	local N=string(e(N),"%9.0fc")
	local R=string(e(r2),"%9.3f")
	local p=2*ttail(e(df_r),abs(_b[`x']/_se[`x']))
	if `p'<=0.1{
		local star "*"
	}
	if `p'<=0.05{
		local  star "**"
	}
	if `p'<=0.01{
		local  star "***"
	}
	if `p'>0.1{
		local star ""
	}


	binscatterhist eduyears `x' [aw=w0] , absorb(w1state) ///
	histogram(`x') xhistbarheight(100) lcolor(gs0) mcolor(edkblue) xcolor(gs6%30) msymbols(O) ymin(13)  ///
	xtitle("`gtitle'", size(medsmall)) ytitle("") title("`title_nc'", size(large) color(gs0)) ///
	text(16.2 -2.2 "Slope: `b'`star' (`se')" "N=`N'" "R{superscript:2}=`R'", place(se) color(gs0) just(left) margin(l+1 t+1 b+1 r+1) fcolor(gs0%0) size(small))
	graph save "${graphpath}corr_nc_`x'.gph", replace


	/*Controlled*/
	reghdfe eduyears `x' $controls [aw=w0], absorb(w1state) cluster(HsId)
	local b=string(_b[`x'],"%9.3f")
	local se=string(_se[`x'],"%9.3f")
	local N=string(e(N),"%9.0fc")
	local R=string(e(r2),"%9.3f")
	local p=2*ttail(e(df_r),abs(_b[`x']/_se[`x']))
	if `p'<=0.1{
		local star "*"
	}
	if `p'<=0.05{
		local  star "**"
	}
	if `p'<=0.01{
		local  star "***"
	}
	if `p'>0.1{
		local star ""
	}


	binscatterhist eduyears `x' [aw=w0] , absorb(w1state) controls($controls) ///
	histogram(`x') xhistbarheight(100) lcolor(gs0) mcolor(edkblue) xcolor(gs6%30)  msymbols(Oh) ymin(13) ///
	xtitle("`gtitle'", size(medsmall)) ytitle("") title("`title_wc'", size(large) color(gs0)) ///
	text(16.2 -2.2 "Slope: `b'`star' (`se')" "N=`N'" "R{superscript:2}=`R'", place(se) color(gs0) just(left) margin(l+1 t+1 b+1 r+1) fcolor(gs0%0) size(small))
	graph save "${graphpath}corr_wc_`x'.gph", replace
}

graph combine 	"${graphpath}corr_nc_pgs_edu" "${graphpath}corr_wc_pgs_edu" ///
				"${graphpath}corr_nc_f_over1" "${graphpath}corr_wc_f_over1"  ///
				, cols(2)  ycommon xcommon l1("Educational Attainment (in Years)") imargin(0 0 0 0)
graph display , xsize(4)
graph export "${graphpath}`name'.`format'", replace

local title "Association of Educational Attainment with \PGS and \IQuali\negspace"
local desc ///
"This figure shows the correlation of completed years of education with" ///
"\PGS and \IQuali\negspace, respectively." ///
"We bin scatterplots using 20 quantiles of the variable of interest." ///
"Gray bars indicate density distributions of the (residualized) variable of interest." ///
"Black lines are fitted from linear regressions of educational attainment on" ///
"the variable of interest." ///
"In the left-column, we only control for state fixed effects." ///
"In the right column, we introduce the full set of control variables."

texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize')  ///
			note_paper("`source' `desc' `controls' `controlfct' `standard' `signif' `error'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})


// -----------------------------------------------------
// Figure 2
// -----------------------------------------------------
local name fig2

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace
xtile f_over1_cat=f_over1, nquantiles(10)
xtile sesbelsky_cat=sesbelsky, nquantiles(10)
keep if inrange(pgs_edu,-4,4)


joyplot pgs_edu, by(f_over1_cat) yrev overlap(10) alpha(50) bwidth(0.5) lc(white%100) lw(0.4) ///
yl ylpattern(dot) ylabsize(3) ylabpos(left) ///
xtitle("PGI{superscript:EA}") ytitle("Decile of Q{superscript:S}")
graph export "${graphpath}`name'.`format'", replace

local title "Distribution of \PGS by \IQuali"
local desc  "This figure shows unconditional \PGS distribution by deciles of" ///
			"\IQuali\negspace. Density distributions are smoothed"  ///
			"using the Epanechnikov kernel function with a bandwidth of 0.5."

texgraph, 	title(`title') name(`name') replace ///
		width(1) file(`format') notesize(`notesize')  ///
		note_paper("`source' `desc'") data(`data')  ///
		inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure 3
// -----------------------------------------------------
local name fig3

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

reghdfe eduyears (c.pgs_edu)##(c.f_over1) ///
$controls [pw=w0], absorb(w1state) cluster(HsId)
margins, at(pgs_edu=(-4.0(.5)3.0) f_over1 =(-2.5(.5)1.5)) saving("${temp}fig_3", replace)

use "${temp}fig_3", clear
colorpalette viridis, nograph
graph twoway ///
(contour _margin _at1 _at2, levels(10) ///
ccolors( "`r(p1)'" "`r(p2)'"  "`r(p3)'"  "`r(p4)'"  "`r(p6)'" ///
		 "`r(p8)'"  "`r(p10)'"  "`r(p12)'"  "`r(p14)'" "`r(p15)'")) ///
, ///
xtitle("Q{superscript:S}", size(large)) xsize(8) ///
ytitle("PGI{superscript:EA}", size(large)) yscale(range(-4(.5)3.0)) ylabel(-4(1)3.0) xscale(range(-2.5(1)1.5)) xlabel(-2.5(1)1.5) ///
ztitle("Education in Years" "(Prediction)") zlabel(, format(%9.1f)) ///
clegend(position(4) bmargin(zero))
graph export "${graphpath}`name'.`format'", replace

local title "Association of Years of Education with \PGS by \IQuali"
local desc  "This figure shows predictions of completed years of education" ///
			"by \PGS and \IQuali cell." ///
			"Predictions are calculated using the model estimated in column (1)" ///
			"of Table \ref{tab: tab3} while allowing for non-linear effects over the range of -4.0(0.5)3.0 of \PGS and -2.5(0.5)1.5 of \IQuali\negspace."

texgraph, 	title(`title')  name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})

// -----------------------------------------------------
// Figure 4
// -----------------------------------------------------
local name fig4

use "${data}data_estim.dta" if sample1==1, clear
qui: gstats transform (standardize)	$s_vars [aw=w0] , replace

postfile fig_4 pos be cu cl str10 ses using "${temp}fig_4.dta", replace
local p=1
foreach y of varlist hs coll2y coll4y postgrad{
	reghdfe `y' (c.pgs_edu)##(c.sesbelsky c.f_over1) ///
	$controls [pw=w0], absorb(w1state) cluster(HsId) keepmata

	foreach r in f_over1 sesbelsky{
		if "`r'"=="f_over1"{
			local pos=`p'-0.15
		}
		if "`r'"=="sesbelsky"{
			local pos=`p'+0.15
		}
		local b=_b[c.pgs_edu#c.`r']
		local cu=_b[c.pgs_edu#c.`r']+_se[c.pgs_edu#c.`r']*invnorm(0.95)
		local cl=_b[c.pgs_edu#c.`r']-_se[c.pgs_edu#c.`r']*invnorm(0.95)

		post fig_4 (`pos') (`b') (`cu') (`cl') ("`r'")

	}
	local ++p
}
postclose fig_4
use "${temp}fig_4.dta", clear

graph twoway ///
(rcap cu cl pos if ses=="f_over1", lcolor("`c1'"%30))  ///
(scatter be pos if ses=="f_over1", mcolor("`c1'") msymbol(O) msize(large))  ///
(rcap cu cl pos if ses=="sesbelsky", lcolor("`c3'"%30))  ///
(scatter be pos if ses=="sesbelsky", mcolor("`c3'") msymbol(T) msize(large)) ///
,  ///
text(-.03 3.5 "School Quality", color("`c1'")) ///
text(+.03 3.5 "Family SES", color("`c3'")) ///
yline(0, lpattern(dash))  ///
ytitle("Coefficient of Child Environment x PGI{superscript:EA}")  ///
xscale(range(1(1)4)) xtitle("") ///
xlabel(1 "High School" 2 "2-year College" 3 "4-year College" 4 "Post-Graduate") ///
legend(off)
graph export "${graphpath}`name'.`format'", replace

local title "Association of Education Degrees with \PGS by \IQuali and Family SES"
local desc  "This figure shows point estimates and 90\% confidence bands for the" ///
			"interaction of \PGS with \IQuali, the interaction of" ///
			"\PGS with an indicator for family SES, and their association with education degrees." ///
			"For each outcome, coefficients are estimated jointly following the specification of equation (\ref{eq: estim_educ})."
texgraph, 	title(`title') name(`name') replace ///
			width(1) file(`format') notesize(`notesize') ///
			note_paper("`source' `desc' `controls' `controlfct' `standard' `error'") data(`data') ///
			inpath(${graphpath}) path_paper(${paper_graphs})


*delete temporary files
local list : dir "${graphpath}"  files "*.gph"
di `list'
foreach f of local list {
    erase "${graphpath}`f'"
}
